%% Skript for setting up parameters for simple physical model
% 
% The state consists of available runoff q, reservoir level V and demand d
% at time t. We use a compact state representation by means of a dimensionless
% number dimL(t+1)=(q(t+1)+V(t+1)/dt)/d(t+1) for dimensionality reduction. 
% h(t) is the reservoir level. The action a(t) is continuous and defined as 
% how much water to consumptively utilize at time t.
% the agent tries to maximize the sum(1/(1+delta)^t*-abs(c-d)) which is defined to be the
% total discounted future (rate delta) benefit. the utility he derives from
% the consumption is denoted by u
% the physical system is such that c is only available from the reservoir.
% reservoir levels are determined by the following balance equation.
% h(t+1) = h(t) + ((q-c)*dt/a)^(1/b);
% any excess water in the reservoir is lost irretrevably.
% 
% Tobias Siegfried, Diego Pontoniero, 16.04.2007
%% Directory
% to change eventually.
% if ispc
%     cd('D:\matlab\diegoGame\intContLearn')
%  elseif ismac
%     cd('/Users/tobias/Documents/Science/matlab/agents/diegoGame/intContLearn/');
% end
%% Parameters
% General
clear all
param.randS = 2; % random seed
param.yearEcho = 1; % plot each xth year on the screen so as to see simulation progress.
% Time
%<<<<<<< .mine
%<<<<<<< .mine
%param.years = 75; % total simulation time
%=======
param.years = 30; % total simulation time
%=======
%param.years = 1; % total simulation time
%>>>>>>> .r105
%>>>>>>> .r96
param.nDays = 365; % number of days per simulation
param.tEnd = param.years * param.nDays; % total umber of timesteps
param.dt = 3600*24; % seconds / decision period.
% Runoff generation (SUPPLY)
param.muS = 0; % runoff generator parameter (mean runoff) [m^3/s]
param.rhoS = .9992; % lag 1 correlation for Markov process
param.sigmaS = 50; % runoff generator parameter (standard deviation to be used for stochastic process)
param.nTS = 1; % number of runoff timeseries to be generated
% Demand generation (DEMAND)
param.muD = 100; % demand generator, assuming muD demand (mean) 
param.sigmaD = 75; % demand generator for seasonal variability
param.dGrowth = .001; % slope of demand growth (linear growth assumed)
% Reservoir parameters
param.aCoeff = .001; % parameter to established reservoir level (h) - volume (V) relationship
param.bCoeff = .5; % ditto
%<<<<<<< .mine
%param.V0 = 10000000000; % Initial volume in m^3
%=======
param.V0 = 0; % Initial volume in m^3
%>>>>>>> .r105
param.Vmax = 80000000000; % Maximum reservoir capacity
% Discount factor
param.delta = .01;
% Learning parameters
param.learnR = .1;
param.decreaseP = 1.0004;
param.explorePIni = .99 ;
%<<<<<<< .mine
param.exploreP = param.explorePIni;
%param.minExpP = 10^-4;
%=======
param.minExpP = 10^-2;
%>>>>>>> .r105
param.discS = 1; % state discretization
%<<<<<<< .mine
param.discQ = 1; % action discretization
param.smoothR = 1;
% Setting up environment
%=======
%param.discQ = 2; % action discretization
% param.maxDimL = 500; % that is the state variable maximum (used for discretization purposes)
% net
param.num_hidden = 40;
param.num_output = 1;
param.min = -1; 
param.max = 1;
% Setting up environment
%>>>>>>> .r105
% Supply
param.q = generateRunoffTS(param);
% Demand
param.d = demandMM(param.muD,param.sigmaD,param.tEnd,param.dGrowth)';
param.scaleF = 5000;
param.dAg = param.scaleF * repmat(diff(tansig(-5:.1/3.65:5)),1,param.years);

%<<<<<<< .mine
% do rest of parameter initialization
%=======
%param.q = q;
%param.d = d;
% do rest of parameter initialization
%>>>>>>> .r105
%param.maxDimL = max(param.q+param.Vmax/param.dt)/min(param.d);
param.logDimL = 0;
if param.logDimL
    param.maxDimL = log(max(param.q+param.Vmax/param.dt)/min(param.d)); % log bsc Doom idea;
    param.minDimL = -10;
else
    param.maxDimL = max(param.q+param.Vmax/param.dt)/min(param.d);
    param.minDimL = 0;
end
param.tBFlow = 0;
%% NO LEARNING
% while assuming that the agent always consumes as much as he can so as to
% % cover the demand.
% V = zeros(param.tEnd+1,1); 
% a = zeros(param.tEnd,1); % action
% l = a;
% r = a; 
% dimL = a;
% V(1) = param.V0;
% for t = 1: param.tEnd
%     if t == 1
%         dimL(t) = abs((q(t) + V(t)/param.dt)) / d(t); % just to get the state at time t=1;
%     end
%     % discuss the availability cases
%     if V(t)/param.dt > d(t)
%         a(t) = d(t);
%     else
%         a(t) = V(t)/param.dt;
%     end
%     V(t+1) = V(t) + q(t)*param.dt - a(t) * param.dt;
%     if V(t+1) < 0,
%         V(t+1) = 0;
%     elseif V(t+1) > param.Vmax
%         l(t) = (V(t+1) - param.Vmax) / param.dt; % l is the loss
%         V(t+1) = param.Vmax;
%     end
%     r(t) = a(t) - d(t);
%     if ~(t==param.tEnd)
%         dimL(t+1) = abs((q(t+1) + V(t+1)/param.dt)) / d(t+1); % new state.
%     end
% end
% h = resVH(V,param.aCoeff,param.bCoeff,'vh');
% V(end) = [];
% % Plot - no learning results
% figure; title('Rule based, no learning')
% subplot(1,2,1);
% plot(q,'b'), hold on; plot(d,'r'); plot(h,'k'); plot(a,'g');hold off;
% xlabel('time');
% subplot(1,2,2)
% semilogy(dimL,'c'); xlabel('time');
% ylabel('state (dimL)');

